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Abstract 

In this paper, a parametric simplex algorithm for solving linear vector optimization prob¬ 
lems (LVOPs) is presented. This algorithm can be seen as a variant of the multi-objective 
simplex (the Evans-Steuer) algorithm [15]. Different from it, the proposed algorithm 
works in the parameter space and does not aim to find the set of all efficient solutions. 
Instead, it finds a solution in the sense of Lohne m, that is, it finds a subset of efficient 
solutions that allows to generate the whole efficient frontier. In that sense, it can also be 
seen as a generalization of the parametric self-dual simplex algorithm, which originally 
is designed for solving single objective linear optimization problems, and is modified to 
solve two objective bounded LVOPs with the positive orthant as the ordering cone in 
Ruszczyhski and Vanderbei m- The algorithm proposed here works for any dimension, 
any solid pointed polyhedral ordering cone C and for bounded as well as unbounded 
problems. 

Numerical results are provided to compare the proposed algorithm with an objective 
space based LVOP algorithm (Benson’s algorithm in [E]), that also provides a solution in 
the sense of [19], and with the Evans-Steuer algorithm m- The results show that for non¬ 
degenerate problems the proposed algorithm outperforms Benson’s algorithm and is on par 
with the Evan-Steuer algorithm. Eor highly degenerate problems Benson’s algorithm m 
outperforms the simplex-type algorithms; however, the parametric simplex algorithm is 
for these problems computationally much more efficient than the Evans-Steuer algorithm. 

Keywords: Linear vector optimization, multiple objective optimization, algorithms, pa¬ 
rameter space segmentation. 

MSC 2010 Classification: 90C29, 90C05, 90-08 


1 Introduction 

Vector optimization problems have been studied for decades and many methods have been 
developed to solve or approximately solve them. In particular, there are a variety of algorithms 
to solve linear vector optimization problems (LVOPs). 
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1.1 Related literature 


Among the algorithms that can solve LVOPs, some are extensions of the simplex method 
and are working in the variable space. In 1973, Evans and Steuer m developed a multi¬ 
objective simplex algorithm that finds the set of all ’efficient extreme solutions’ and the set 
of all ’unbounded efficient edges’ in the variable space, see also [HI Algorithm 7.1]. Later, 
some variants of this algorithm have been developed, see for instance [H El El uni CZl E3]. 
More recently, Ehrgott, Puerto and Rodriguez-Chia m developed a primal-dual simplex 
method that works in the parameter space. This algorithm does not guarantee to find the set 
of all efficient solutions, instead it provides a subset of efficient solutions that are enough to 
generate the whole efficient frontier in case the problem is ’bounded’. All of these simplex-type 
algorithms are designed to solve LVOPs with any number of objective functions where the 
ordering is component-wise. Among these, the Evans-Steuer algorithm m is implemented 
as a software called ADBASE m- The idea of decomposing the parameter set is also used 
to solve multiobjective integer programs, see for instance [26] . 

In [27], Ruszczyhski and Vanderbei developed an algorithm to solve LVOPs with two 
objectives and the efficiency of this algorithm is equivalent to solving a single scalar linear 
program by the parametric simplex algorithm. Indeed, the algorithm is a modification of the 
parametric simplex method and it produces a subset of efficient solutions that generate the 
whole efficient frontier in case the problem is bounded. 

Apart from the algorithms that work in the variable or parameter space, there are al¬ 
gorithms working in the objective space. In [8], Dauer and Liu proposed a procedure to 
determine the ’maximal’ extreme points and edges of the image of the feasible region. Later, 
Benson [5] proposed an outer approximation algorithm that also works in the objective space. 
These methods are motivated by the observation that the dimension of the objective space is 
usually much smaller than the dimension of the variable space, and decision makers tend to 
choose a solution based on objective values rather than variable values, see for instance [7]. 
Lohne m introduced a solution concept for LVOPs that takes into account these ideas. Ac¬ 
cordingly a solution consists of a set of ’point maximizers (efficient solutions)’ and a set of 
’direction maximizers (unbounded efficient edges)’, which altogether generate the whole effi¬ 
cient frontier. If a problem is ’unbounded’, then a solution needs to have a nonempty set of 
direction maximizers. There are several variants of Benson’s algorithm for LVOPs. Some of 
them can also solve unbounded problems as long as the image has at least one vertex, but only 
by using an additional Phase 1 algorithm, see for instance m Section 5.4]. The algorithms 
provided in na nni ESI EDI solve in each iteration at least two LPs that are of the same size 
as the original problem. An improved variant where only one LP has to be solved in each 
iteration has been proposed independently in [6] and m- In addition to solving (at least) 
one LP, these algorithms solve a vertex enumeration problem in each iteration. As it can be 
seen in iniE2lE3ES!> it is also possible to employ an online vertex enumeration method. In 
this case, instead of solving a vertex enumeration problem from scratch in each iteration, the 
vertices are updated after an addition of a new inequality. Recently, Benson’s algorithm was 
extended to approximately solve bounded convex vector optimization problems in [HEI]. 

1.2 The proposed algorithm 

In this paper, we develop a parametric simplex algorithm to solve LVOPs of any size and 
with any solid pointed polyhedral ordering cone. Although the structure of the algorithm is 
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similar to the Evans-Steuer algorithm, it is different since the algorithm proposed here works 
in the parameter space and it finds a solution in the sense that Lohne proposed in [19]. In 
other words, instead of generating the set of all point and direction maximizers, it only finds a 
subset of them that already allows to generate the whole efficient frontier. More specifically, 
the difference can be seen at two points. First, in each iteration instead of performing a 
pivot for each ’efficient nonbasic variable’, we perform a pivot only for a subset of them. 
This already decreases the total number of pivots performed throughout the algorithm. In 
addition, the method of finding this subset of efficient nonbasic variables is more efficient 
than the method that is needed to find the whole set. Secondly, for an entering variable, 
instead of performing all possible pivots for all ’efficient basic variables’ as in |15| . we perform 
a single pivot by picking only one of them as the leaving variable. In this sense, the algorithm 
provided here can also be seen as a generalization of the algorithm proposed by Ruszczyhski 
and Vanderbei m to unbounded LVOPs with more than two objectives and with more 
general ordering cones. 

In each iteration the algorithm provides a set of parameters which make the current vertex 
optimal. This parameter set is given by a set of inequalities among which the redundant ones 
are eliminated. This is an easier procedure than the vertex enumeration problem, which is 
required in some objective space algorithms. Different from the objective space algorithms, 
the algorithm provided here does not require to solve an additional LP in each iteration. 
Moreover, the parametric simplex algorithm works also for unbounded problems even if the 
image has no vertices and generates direction maximizers at no additional cost. 

As in the scalar case, the efficiency of simplex-type algorithms is expected to be better 
whenever the problem is non-degenerate. In vector optimization problems, one may observe 
different types of redundancies if the problem is degenerate. The first one corresponds to the 
’primal degeneracy’ concept in scalar problems. In this case, a simplex-type algorithm may 
find the same ’efficient solution’ for many iterations. That is to say, one remains at the same 
vertex of the feasible region for more than one iteration. The second type of redundancy 
corresponds to the ’dual degeneracy’ concept in scalar problems. Accordingly, the algorithm 
may find different ’efficient solutions’ which yield the same objective values. In other words, 
one remains at the same vertex of the image of the feasible region. Additionally to these, a 
simplex-type algorithm for LVOPs may find efficient solutions which yield objective values 
that are not vertices of the image of the feasible region. Note that these points that are on 
a non-vertex face of the image set are not necessary to generate the whole efficient frontier. 
Thus, one can consider these solutions also as redundant. 

The parametric simplex algorithm provided here may also find redundant solutions. How¬ 
ever, it will be shown that the algorithm terminates at a finite time, that is, there is no 
risk of cycling. Moreover, compared to the Evans-Steuer algorithm, the parametric simplex 
algorithm finds much fewer redundant solutions in general. 

We provide different initialization methods. One of the methods requires to solve two 
LPs while a second method can be seen as a Phase 1 algorithm. Both of these methods work 
for any LVOP. Depending on the structure of the problem, it is also possible to initialize the 
algorithm without solving an LP or performing a Phase 1 algorithm. 

This paper is structured as follows. Sectionis dedicated to basic concepts and notation. 
In Sectionj^ the linear vector optimization problem and solution concepts are introduced. The 
parametric simplex algorithm is provided in Section]^ Different methods of initialization are 
explained in Section [T8| Illustrative examples are given in Sectionj^ In Sectionj^ we compare 
the parametric simplex algorithm provided here with the different simplex algorithms for 
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solving LVOPs that are available in the literature. Finally, some numerical results regarding 
the efficiency of the proposed algorithm compared to Benson’s algorithm and the Evans-Steuer 
algorithm are provided in Section 

2 Preliminaries 

For a set ^ C , int A, riA, cl A, bdA, convA, cone A denote the complement, interior, 

relative interior, closure, boundary, convex hull, and conic hull of it, respectively. If A C 
is a non-empty polyhedral convex set, it can be represented as 

A = conv ..., x^} + cone {A:^,..., (1) 

where s G N \ {0},t G N, each G is a point, and each y G \ {0} is a direction of 
A. Note that A: G \ {0} is called a direction of A if A -h {ak G M^| a > 0} C A. The set 
Aoo := cone {A:^,..., A:^} is the recession cone of A. The set of points ..., x^} together 
with the set of directions {A:^,..., A:^} are called the generators of the polyhedral convex set 
A. We say ({x^,..., {A:^,..., A;^}) is a V-representation of A whenever 0 holds. For 

convenience, we define cone0 = {0}. 

A convex cone C is said to be non-trivial if {0} C (7 C and pointed if it does not 
contain any line. A non-trivial convex pointed cone C defines a partial ordering <c on 

V <c w w — V ^ C. 

For a non-trivial convex pointed cone C C a point y ^ A is called a C-maximal element 
of A if {{y} + C \ {0}) A — tj). If the cone C is solid^ that is, if it has a non-empty 
interior, then a point ^ G A is called weakly C-maximal if {{y} -h int C) H A = 0. The set 
of all (weakly) (7-maximal elements of A is denoted by (w)Max^(A). The set of (weakly) 
C-minimal elements is defined by (w)Min^ (A) := (w)Max_^ (A). The (positive) dual cone 
of C is the set [z ^W^\\/y ^ C : z^y > O}. The positive orthant of is denoted by 

M^, that is, := {y G M^| > 0, z = 1,..., g}. 

3 Linear vector optimization problems 

We consider a linear vector optimization problem (LVOP) in the following form 

maximize P^x (with respect to <c) (P) 

subject to Ax < 6, 

X > 0, 

where P G A G b G M^, and (7 C is a solid polyhedral pointed ordering 

cone. We denote the feasible set by A := {x G W^\ Ax <b^x> 0}. Throughout, we assume 
that 0 is feasible, i.e., A ^ 0. The image of the feasible set is defined as P^[A] := {P^x G 
M^l X G A}. 

We consider the solution concept for LVOPs as in m- To do so, let us recall the following. 
A point x G A is said to be a (weak) maximizer for if X is (weakly) C-maximal in 
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P[A’]. The set of (weak) maximizers of (|^ is denoted by (w)Max([^. The homogeneous 
problem of 0 is given by 

maximize x (with respect to <c) (P^) 

subject to Ax < 0, 

X > 0. 


The feasible region of (P^), namely {x G W^\ Ax < 0, x > 0}, satisfies X^ = Tbo, that 

is, the non-zero points in X^ are exactly the directions of X. A direction k G \ {0} of 
X is called a (weak) maximizer ioY 0 i f th e corresponding point k G X^ \ {0} is a (weak) 
maximizer of the homogeneous problem (P^). 


Definition 3.1 ([lEKIH]). A set X C X is ealled a set of feasible points for 0 and a set 
X^ C X^ \ {0} is ealled a set of feasible directions for 0. 

A pair of sets (T, T^) is ealled a finite supremizer for 0 ifX is a non-empty finite set 
of feasible points for 0. is a (not neeessarily non-empty) finite set of feasible direetions 
for Q, and 

conv P^[X] + cone P^[X^] -C = P^[X]- C. (2) 

A finite supremizer (X, X^) of 0 is ealled a solution to 0 if it eonsists of only maxi¬ 
mizers. 


The set V := P^[T] — C is called the lower image of Let , 2 /^ be the generating 

vectors of the ordering cone C. Then, ({0} , ..., y^}) is a V-representation of the cone (7, 

that is, C = cone ..., y^}. Clearly, if (A, X^) is a finite supremizer, then (P^[T], 

..., —y^}) is a V-representation of the lower image V. 

Definition 3.2. Q is said to be bounded if there exists p G sueh that V C {p} — C. 

Remark 3.3. Note that the recession cone of the lower image, is equal to the lower 
image of the homogeneous problem, that is, V^o — — C, see m Lemma 4.61]. Clearly, 

Poo 3 —C, which also implies P^ C —C^. In particular, if Q is bounded, then we have 
p^ = -C and X^ = 0. 

The weighted sum scalarized problem for a parameter vector w G is 

maximize w^P^x (PiH) 

subject to Ax < b, 

X >0, 


and the following well known proposition holds. 


Proposition 3.4 ([211 Theorem 2.5]). A point x ^ X is a maximizer (weak maximizer) 
0/0 if and only if it is an optimal solution to ( |Pi(re)[ ) for some w G int (u.eC+\{0}). 

Proposition |3.4| suggests that if one could generate optimal solutions, whenever they exist, 
to the problems (Pi (re) ) for re G int then this set of optimal solutions X would be a set of 
(point) maximizers of 0. Indeed, it will be enough to solve problem (EMI for re G ri VP, 
where 


VP := {re G C^\ w^c = 1}, 


(3) 
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for some fixed c G int C. Note that ( |Pi(i(;)D is not necessarily bounded for all w ^ li W. 
Denote the set of all w ^ liW such that ( |Pi(r(;)| ) has an optimal solution by If one can 
find a finite partition of VI 4 such that for each i G {1,... ,s} there exists x'^ ^ X 

which is an optimal solution to ( |Pi(r(;)[ ) for all w G then, clearly, X — ..., will 

satisfy (§ provided one can also generate a finite set of (direction) maximizers X^. Trivially, 
if problem 0 is bounded, then ( |Pi(re)[ ) can be solved optimally for all w G X^ — 0, 
and (T, 0) satisfies (§. If problem (|P[) is unbounded, we will construct in Section]^ a set X^ 
by adding certain directions to it whenever one encounters a set of weight vectors w G 
for which ( Pi(r(;)| ) cannot be solved optimally. The following proposition will be used to 
prove that this set together with X — will indeed satisfy (§. It provides 

a characterization of the recession cone of the lower image in terms of the weighted sum 
scalarized problems. More precisely, the negative of the dual of the recession cone of the 
lower image can be shown to consist of those w G for which ( |Pi(r(;) ) can be optimally 
solved. 

Proposition 3.5. The recession cone Voo of the lower image satisfies 

—V^ = {le G C^\ ( |Pi(r(;)D is bounded}. 
we have Voo = — C. Using 0 G X^ and 0 G C, we obtain 


Proof By Remark 
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-v+ ={we M'^l eX^,yceC: vF{P'^x^ - c) < 0} 

^{weC^l : vFp'^x^ < 0 }. 


Let w G —Vtoi consider the weighted sum scalarized problem of (P^) given by 


maximize 
subject to 


P^x 
Ax < 0, 

X > 0. 


(4) 


(F'liw)) 


By®, = 0 is an optimal solution, which implies by strong duality of the linear program 
that there exist G with > Pw and y* > 0. Then, is also dual feasible 

for ( |Pi(ie)| ). By the weak duality theorem, ( |Pi(r(;)[ ) can not be unbounded. 

For the reverse inclusion, let w G be such that ( |Pi(re)| ) is bounded, or equivalently, 
an optimal solution exists for ( |Pi(ie) ) as we assume T ^ 0. By strong dual ity, the dual 
problem of ( Pi(r(;)[) h as an optimal solution ^*, which is also dual feasible for (Pi{w)). By 
weak duality, (Pi {w)) is bounded and has an optimal solution x^. Then, w G —V^ holds. 
Indeed, assuming the contrary, one can easily find a contradiction to the optimality of □ 


4 The parametric simplex algorithm for LVOPs 


In [15], Evans and Steuer proposed a simplex algorithm to solve linear multiobjective opti¬ 
mization problems. The algorithm moves from one vertex of the feasible region to another 
until it finds the set of all extreme (point and direction) maximizers. In this paper we pro¬ 
pose a parametric simplex algorithm to solve LVOPs where the structure of the algorithm is 
similar to the Evans-Steuer algorithm. Different from it, the parametric simplex algorithm 


provides a solution in the sense of Definition 3.1, that is, it finds subsets of extreme point and 


6 



































direction maximizers that generate the lower image. This allows the algorithm to deal with 
the degenerate problems more efficiently than the Evans-Steuer algorithm. More detailed 
comparison of the two algorithms can be seen in Section 

In [27], Ruszczyhski and Vanderbei generalize the parametric self dual method, which 
originally is designed to solve scalar LPs [32], to solve two-objective bounded linear vector 
optimization problems. This is done by treating the second objective function as the auxiliary 
function of the parametric self dual algorithm. The algorithm provided here can be seen as 
a generalization of the parametric simplex algorithm from biobjective bounded LVOPs to q- 
objective LVOPs {q > 2) that are not necessarily bounded where we also allow for an arbitrary 
solid polyhedral pointed ordering cone C. 

We first explain the algorithm for problems that have a solution. One can keep in mind 


that the methods of initialization proposed in Section 4.8 will verify if the problem has a 
solution or not. 

Assumption 4.1. There exists a solution to problem 

This assumption is equivalent to having a nontrivial lower image 'P, that is, 0 ^ P C 
Clearly P ^ 0 implies T ^ 0, which is equivalent to our standing assumption. Moreover, 
by Definition |3.1| and Proposition |3.4| Assumption |4.1| implies that there exists a maximizer 
which guarantees that there exists some G intC^ such that problem has an 

optimal solution. In Section 4^, we will propose methods to find such a It will be seen 
that the algorithm provided here finds a solution if there exists one. 


4.1 The parameter set A 

Throughout the algorithm we consider the scalarized problem (EMI' for w ^ W where W 
is given by (© for some fixed c G int (7. As VP is g — 1 dimensional, we will transform the 
parameter set W into a set A C Assume without loss of generality that Cq = 1. Indeed, 

since C was assumed to be a solid cone, there exists some c G int (7 such that either Cq — I 
OT Cq = —1. For Cg = — 1, one can consider problem 0 where (7 and P are replaced by —(7 
and —P. 

Let c = (ci,..., Cg_i)^ G and define the function w{X) : ^ and the set 

A C as follows: 

^(A) := (Ai,..., A^_i, 1 - c^A)^, 

A := {A G w{X) G (7+}. 

As we assume Cq = 1, c^w{X) = 1 holds for all A G A. Then, w{X) G W for all A G A 
and for any w G VP, {wi,... ,Wq-i)^ G A. Moreover, if A G intA, then w{X) G riVP and 
a w ^ riVP, then {wi^... ^Wq-i)^ G intA. Throughout the algorithm, we consider the 
parametrized problem 

(Pa) := (Pi(^A))) 

for some generic A G 


4.2 Segmentation of A: dictionaries and their optimality region 

We will use the terminology for the simplex algorithm as it is used in [32]. First, we introduce 
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slack variables [x^+i,..., Xn+rnV obtain x G and rewrite (Pa) as follows 

maximize w{X)^[P^ 0]x (Pa) 

subject to [A I]x = 6, 

X > 0, 

where / is the identity and 0 is the zero matrix, all in the correct sizes. We consider a partition 
of the variable indices {1, 2,..., n + m} into two sets B and J\f. Variables xj^j G B^ are called 
basic variables and xj^j G J\f^ are called nonbasic variables. We write x = 
permute the columns of [A I] to obtain [B V] satisfying [A I]x — Bxs + Nxj\f^ where 
B G and N G Similarly, we form matrices Ps G and Pj\f G such 

that [P^0]x = P^XB + PjfX^. In order to keep the notation simple, instead of writing [P^ 0]x 
we will occasionally write where x stands then for the original decision variables in 

without the slack variables. 

Whenever B is nonsingular, xs can be written in terms of the nonbasic variables as 
xs — B~^b — B~'^Nxj^. Then, the objective function of (Pa) is 

w{X)^[P^ 0]x = ic(A)^^(A) — w{XY'ZjfX^.^ 

where ^(A) = P^B-^b and = {B-^N)^Pb - PaT- 

We say that each choice of basic and nonbasic variables defines a unique dictionary. 
Denote the dictionary defined by B and J\f hy D. The basic solution that corresponds to D is 
obtained by setting xj\f = 0. In this case, the values of the basic variables become B~^b. Both 
the dictionary and the basic solution corresponding to this dictionary are said to be primal 
feasible if B~^b > 0. Moreover, if w{X)^Zjf > 0, then we say that the dictionary D and the 
corresponding basic solution are dual feasible. We call a dictionary and the corresponding 
basic solution optimal if they are both primal and dual feasible. 

For j G A/", introduce the halfspace 

If {A € R^-^\ w{Xfzfe^ > 0}, 

G denotes the unit column vector with the entry corresponding to the variable 
1. Note that if D is known to be primal feasible, then D is optimal for A G A^, 

A® := fl If. 

jeJ\f 

The set H A is said to be the optimality region of dictionary D. 

Proposition |3.4| already shows that a basic solution corresponding to a dictionary D with 
A^ n A ^ 0 yields a (weak) maximizer of 0 Throughout the algorithm we will move from 
dictionary to dictionary and collect their basic solutions into a set X. We will later show that 
this set will be part of the solution (T, X^) of 0. The algorithm will yield a partition of the 
parameter set A into optimality regions of dictionaries and regions where (Pa) is unbounded. 
The next subsections explain how to move from one dictionary to another and how to detect 
and deal with unbounded problems. 


where 
Xj being 
where 


4.3 The set of entering variables 

We call a defining (non-redundant) collection of half-spaces of the optimality region 

n A if C J\f satisfies 

A^ n A = n A and 

(5) 

n A c Pi n A, for any J C . 
j&J 

For a dictionary any nonbasic variable xj^j G J^, is a candidate entering variable. Let us 
call the set an index set of entering variables for dictionary D. 

For each dictionary throughout the algorithm, an index set of entering variables is found. 
This can be done e.g. by the following two methods. Firstly, using the duality of polytopes, the 
problem of finding defining inequalities can be transformed to the problem of finding a convex 
hull of given points. Then, the algorithms developed for this matter, see for instance [3], can 
be employed. Secondly, in order to check if j G A/" corresponds to a defining or a redundant 
inequality one can consider the following linear program in A G 

maximize w{\Y' Zjfe^ 

subject to w{\Y'Zjfe^ > 0, for all j ^ J\f\ {{j} U 
w{XfY > 0, 

where is the index set of redundant inequalities that have been already found and 

Y — ... is the matrix where y^^... ^y^ are the generating vectors of the ordering cone 

C. The inequality corresponding to the nonbasic variable Xj is redundant if and only if an 
optimal solution to this problem yields w{X'Y'Zjj-e^ < 0. In this case, we add j to the set 
jredun^ Otherwise, it is a defining inequality for the region A^ H A and we add j to the set 
. The set is obtained by solving this linear program successively for each untested 
inequality against the remaining. 

Remark 4.2. For the numerical examples provided in Section the second method is em¬ 
ployed. Note that the number of variables for each linear program is g — 1, which is much 
smaller than the number of variables n of the original problem in general. Therefore, each 
linear program can be solved accurately and fast. Thus, this is a reliable and sufficiently 
efficient method to find . 

Before applying one of these methods, one can also employ a modified Fourier-Motzkin 
elimination algorithm as described in |1] in order to decrease the number of redundant in¬ 
equalities. Note that this algorithm has a worst-case complexity of — l)^)n^). Even 

though it does not guarantee to detect all of the redundant inequalities, it decreases the 
number significantly. 

Note that different methods may yield a different collection of indices as the set might 
not be uniquely defined. However, the proposed algorithm works with any choice of . 

4.4 Pivoting 

In order to initialize the algorithm one needs to find a dictionary for the parametrized 
problem (P;^) such that the optimality region of satisfies A^ n int A ^ 0. Note that the 
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existence of is guaranteed by Assumption 4.1 and by Proposition 3.4 There are different 
methods to find an initial dictionary and these will be discussed in Section |4.8[ For now, 
assume that is given. By Proposition 3^, the basic solution corresponding to is a 
maximizer to (^. As part of the initialization, we find an index set of entering variables 
as defined by (^. 

Throughout the algorithm, for each dictionary D with given basic variables optimality 
region A^nA, and index set of entering variables J^, we select an entering variable Xj^j G J^, 
and pick analog to the standard simplex method a leaving variable xi satisfying 


zG argmm . ^ . : 

ieB \B JM)ij 
> 0 


( 6 ) 


whenever there exists some i with > 0. Here, indices i, j are written on behalf of the 

entries that correspond to the basic variable Xi and the nonbasic variable Xj, respectively. Note 
that this rule of picking leaving variables, together with the initialization of the algorithm with 
a primal feasible dictionary guarantees that each dictionary throughout the algorithm is 
primal feasible. 

If there exists a basic variable Xi with > 0 satisfying we perform the pivot 

Xj GG Xi to form the dictionary D with basic variables B — and nonbasic variables 

M — (A/’U{i})\{j}. For dictionary we have Ip — cl [if)^ = {A G w{\Y' Zpe^ < 0}. 

If dictionary D is considered at some point in the algorithm, it is known that the pivot 
Xi GG Xj will yield the dictionary D considered above. Thus, we call (i, j) an explored pivot 
(or direetion) for D. We denote the set of all explored pivots of dictionary D by . 


4.5 Detecting unbounded problems and constructing the set 

Now, consider the case where there is no candidate leaving variable for an entering variable 
Xj, j G J^, of dictionary that is, {B~^Nep < 0. As one can not perform a pivot, it is not 
possible to go beyond the halfspace ip. Indeed, the parametrized problem (P^) is unbounded 
for X ^ ip. The following proposition shows that in that case, a direction of the recession 


cone of the lower image can be found from the current dictionary T), see Remark 3.3 


Proposition 4.3. Let D he a dietionary with basie and nonbasie variables B and J\f, H A 
be its optimality region satisfying A^ Hint A 0, and he an index set of entering variables. 
If for some j G , {B~^Nep < 0, then the direetion x^ defined by setting x^ = —B~^Ne^ 
and x^ = e^ is a maximizer to 0 


and P^x^ = -Zjj-eP 


Proof. Assume {B ^Nep < 0 for j G and define x^ by setting x^ = —B ^Ne^ and 
x^ — ' By definition, the direction x^ would b e a maximizer for 0 if and only if it is 

(point) maximizer for the homogeneous problem (P^), see section It holds 


IS a 


[yl l]x^ = Bx% + Nx^^ = 0. 


Moreover, = e^ > 0 and x^ = —B ^Ne^ > 0 by assumption. Thus, x^ is primal feasible 


for problem (P^) and also for problem (P^(ie(A))) for all A G A, that is, x^ G \ {0}. Let 


A G A^ Hint A, which implies w{X) G ri VF C int C^. Note that by definition of the optimality 
region, it is true that w{X)^Zp > 0. Thus, x^ is also dual feasible for (Pi{w{X))) and it 
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is an optim al so lution for the parametrized homogeneous problem for A G n int A. By 
Proposition 3.4 (applied to (P^) and is a maximizer of (P^). The value of 


the objective function of (P^) at is given by 


[P^' 0]x'^ = = -Zlfe 




□ 


Remark 4.4. If for an entering variable Xj,j G J^, of dictionary D, there is no candidate 
leaving variable, we conclude that problem 0 is unbounded in the sense of Definition |3.2[ 
Then, in addition to the set of point maximizers T one also needs to find the set of (direction) 
maximizers of 0. which by Proposition [4.3| can be obtained by collecting directions x^ 


defined by x^ = and x^ = for every j G with < 0 for all dictionaries 

visited throughout the algorithm. For an index set of entering variables of each dictionary 
D, we denote the set of indices of entering variables with no candidate leaving variable for 
dictionary D by := {j G J^\ B~^Ne^ < 0}. In other words, C is such that for 
any j G J^, (Pa) is unbounded for A ^ 


4.6 Partition of A: putting it all together 

We have seen in the last subsections that basic solutions of dictionaries visited by the al¬ 
gorithm yield (weak) point maximizers of 0 and partition A into optimality regions for 
bounded problems (Pa), while encountering an entering variable with no leaving variable in 
a dictionary yields direction maximizers of © as well as regions of A corresponding to un¬ 
bounded problems (Pa)- This will be the basic idea to construct the two sets X and X^ 
and to obtain a partition of the parameter set A. In order to show that (T, X^) produces 
a solution to 0. one still needs to ensure finiteness of the procedure, that the whole set 
A is covered, and that the basic solutions of dictionaries visited yield not only weak point 
maximizers of 0. but point maximizers. 

Observe that whenever Xj, j G J^, is the entering variable for dictionary D with A^nA ^ 
0 and there exists a leaving variable the optimality region A^ for dictionary D after the 
pivot is guaranteed to be non-empty. Indeed, it is easy to show that 

0 C A^ n A^ C {A G w{XfZxe^ = 0}, 

where J\f is the collection of nonbasic variables of dictionary D. Moreover, the basic solutions 
read from dictionaries D and D are both optimal solutions to the parametrized problem (Pa) 
for A G A^ n A^. Note that the common optimality region of the two dictionaries has no 
interior. 

Remark 4.5. a. In some cases it is possible that A^ itself has no interior and it is a subset 
of the neighboring optimality regions corresponding to some other dictionaries. 

b. Even though it is possible to come across dictionaries with optimality regions having empty 
interior, for any dictionary D found during the algorithm A^ H int A ^ 0 holds. This is 
guaranteed by starting with a dictionary satisfying A^° H int A ^ 0 together with 
the rule of selecting the entering variables, see (§. More specifically, throughout the 
algorithm, whenever corresponds to the boundary of A it is guaranteed that j ^ . 

By this observation and by Proposition [3^ it is clear that the basic solution corresponding 
to the dictionary D is not only a weak maximizer but it is a maximizer. 


11 












Let us denote the set of all primal feasible dictionaries D satisfying H int A ^ 0 by V. 
Note that P is a finite collection. Let the set of parameters A G A yielding bounded scalar 
problems (Pa) be A^. Then it can easily be shown that 


A 5 := {A G A| (Pa) has an optimal solution} (7) 

Dev 

Note that not all dictionaries in V are required to be known in order to cover A^. First, 
the dictionaries mentioned in Remark |4.5| a. do not provide a new region within A^. One 
should keep in mind that the algorithm may still need to visit some of these dictionaries 
in order to go beyond the optimality region of the current one. Secondly, in case there are 
multiple possible leaving variables for the same entering variable, instead of performing all 
possible pivots, it is enough to pick one leaving variable and continue with this choice. Indeed, 
choosing different leaving variables leads to different partitions of the same subregion within 
Kb. 

By this observation, it is clear that there is a subcollection of dictionaries V C V which 
defines a partition of A^ in the following sense 


\J{A^nA) = A,. ( 8 ) 

Dev 


If there is at least one dictionary D with ^ 0, it is 
unbounded. If further A^ is connected, one can show that 


known by Remark 


4.4 


that 


0 is 


f] (ifnA) = Afe, 
Def>,jeJ^ 


holds. Indeed, connectedness of A^ is correct, see Remark 4.7 below. 


(9) 


4.7 The algorithm 

The aim of the parametrized simplex algorithm is to visit a set of dictionaries T) satisfying Q . 
In order to explain the algorithm we introduce the following definition. 


Definition 4.6. D ^ V is said to he a boundary dictionary if hP and an index set of 
entering variables is known. A boundary dietionary is said to be visited if the resulting 
dietionaries of a ll po ssible pivots from D are boundary and the index set eorresponding to 
(see Remark 4-4) known. 


The motivation behind this definition is to treat the dictionaries as nodes and the possible 
pivots between dictionaries as the edges of a graph. Note that more than one dictionary may 
correspond to the same maximizer. 


Remark 4.7. The graph described above is not necessarily connected. However, there exists 
a connected subgraph which includes at least one dictionary corresponding to each maximizer 
found by visiting the whole graph. The proof for the case C = is given in |28j and it can 
be generalized easily to any polyhedral ordering cone. Note that this implies that the set A^ 
is connected. 
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The idea behind the algorithm is to visit a sufficient subset of ’nodes’ to cover the set 
A^. This can be seen as a special online traveling salesman problem. Indeed, we employ the 
terminology used in m- The set of all ’currently’ boundary and visited dictionaries through 
the algorithm are denoted by BD and VS^ respectively. 

The algorithm starts with BD — {D^} and VS — where is the initial dictionary 
with index set of entering variables J^°. We initialize as the empty set and X as 
where is the basic solution corresponding to . Also, as there are no explored directions 
for we set = 0. 

For a boundary dictionary we consider each j G and check the leaving variable 


corresponding to x 


j • 


If there is no leaving vari able, we add x^ defined by x^ — —B ^ 

Otherwise, a corresponding leaving variable 


4.3 


and x^ = to the set see Proposition 
Xi is found. If (j,i) ^ ^ we perform the pivot Xj CG Xi as it has not been explored before. 

We check if the resulting dictionary D is marked as visited or boundary, li D there 

is no need to consider D further. If 5 G BD^ then (i, j) is added to the set of explored 
directions for D. In both cases, we continue by checking some other entering variable of D. 
If D is neither visited nor boundary, then the corresponding basic solution x is added to the 
set T, an index set of entering variables is computed, (i, j) is added to the set of explored 
directions E^ ^ and D itself is added to the set of boundary dictionaries. Whenever all j G 
have been considered, D becomes visited. Thus, D is deleted from the set BD and added to 
the set VS. The algorithm stops when there are no more boundary dictionaries. 

Theorem 4.8. Algorithm^returns a solution (T,T^) to 0. 

Proof. Algorithm [T] terminates in a finite number of iterations since the overall number of 
dictionaries is finite and there is no risk of cycling as the algorithm never performs ’already 
explored pivots’, see line 13. X^X^ are finite sets of feasible points and directions, respec¬ 
tively, f or d^ , and they consist of only maximizers by Proposition s |3.4| and 4.3 together with 
Remark [ 4 . 5| b. Hence, it is enough to show that satisfies Q. 

Observe that by construction, the set of all visited dictionaries V VS Sit termination 
satisfies ([s]). Indeed, there are finitely many dictionaries and A^ is a connected set, see 
Remark 4.7 It is guaranteed by §)_ that for any w G (7^, for which ( |Pi(r(;)| ) is bounded. 


there exists an optimal solution x G T of (Pi{w)). Then, it is clear that (T, X^) satisfies (§ 
as long as R := coneP^[T^] — C is the recession cone Voo of the lower image. 

If for all D the set = 0, then Q is bounded, X^ = 0, and trivially R = —C = Voo- 
For the general case, we show that —V^ = —R^ which implies Voo — cone — C. 


4.4 


and 


4.7 


we 


Assume there is at least one dictionary D with ^ 0. Then, by Remarks 
0 is unbounded, X^ ^ 0 and V also satisfies Q. On the one hand, by definition of ^ 
can write © as 

Afo = P {A G A| w{\f'> 0}, (10) 

D&VS, 

where JVd is the set of nonbasic variables corresponding to dictionary D. On the other hand, 
by construction and by Proposition |4.3| we have 


R = cone {{-j 


u 

Devs 


p} u {- 


-y 


V}), 


where {y ^,..., y*} is the set of generating vectors for the ordering cone C. The dual cone can 


13 












Algorithm 1 Parametric Simplex Algorithm for LVOP 


1 : Find and an index set of entering variables J^°; 
BD = {L>0},A = {a;°}; 


2 : Initialize 


?D° 


,R = 


VS, 

while BD ^ 0 do 

Let D ^ BD with nonbasic variables Af and index set of entering variables J 

for j G do 

Let Xj be the entering variable; 
if B~^Ne^ < 0 then 


Let be such that x^ = —B and x^ = e 


J- 


9 

10 

11 

12 

13 

14 

15 

16 

17 

18 

19 

20 
21 
22 

23 

24 

25 

26 

27 

28 

29 

30 

31 

32 


P^[X^] ^ P^[X^] U 

else ^ 

Pick i € argmin^ge^ {B-^N)ij>o 
if (j, i) ^ then 

Perform the pivot with entering variable xj and leaving variable xf, 
Call the new dictionary D with nonbasic variables M = Af U {i} \ {j}; 
D ^VS then 
if D ^ BD then 

else 

Let X be the basic solution for D] 

X i — X U 

p^[x]^p^[x]yj{p^xy, 

Compute an index set of entering variables of D; 

Let E^^{{i,j)y, 

BD ^ BDU{D}; 

end if 
end if 
end if 
end if 
end for 

VS ^VSU {D}, BD^BD\ {D}; 

end while 


D. 


33 


return 


^X X^^ 


: A finite solution of ©; 
: V representation of V. 
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be written as 


k 

R+= Pi {w eR‘^\w'^Z]^^e^ <0}nf]{w <0}. 

Devs, jeJ^ *=i 


( 11 ) 


3.5 


( Pi(r(;)| ) has an optimal solution. As c^w > 0, 
has an optimal solution, where A := ... ^Wq-i)^ and thus 

G A^. Then, by (10), A G A and 
w{\Y' > 0 for all j G J^, D ^VS. This holds if and only if w ^ ~R+ by definition of 
A and by The other inclusion can be shown symmetrically. □ 


Now, let w G —V^. By proposition 
also (Pi(^)) = (Pa 
w{X) — By the definition of A^ given by ([^, A 


Remark 4.9. In general, simplex-type algorithms are known to work better if the problem 
is not degenerate. If the problem is degenerate. Algorithm may find redundant maximizers. 
The effects of degeneracy will be provided in more detail in SectionFor now, let us mention 
that it is possible to eliminate the redundancies by additional steps in Algorithm [T] There 
are two types of redundant maximizers. 

a. Algorithm may find multiple point (direction) maximizers that are mapped to the same 
point (direction) in the image space. In order to find a solution that is free of these type of 
redundant maximizers, one may change line 21 (9) of the algorithm such that the current 
maximizer x (x^) is added to the set X (X^) only if its image is not in the current set 
P^[X] {P^[Xh]). 


b. Algorithm [T] may also find maximizers whose image is not a vertex on the lower image. 
One can eliminate these maximizers from the set X by performing a vertex elimination at 
the end. 


4.8 Initialization 

There are different ways to initialize Algorithm We provide two methods, both of which 
also determine if the problem has no solution. Note that 0 has no solution if A = 0 or if the 
lower image is equal to that is, if ( |Pi(re)| ) is unbounded for all w G int C^. We assume X 
is nonempty. Moreover, for the purpose of this section, we assume without loss of generality 
that b > 0. Indeed, if 6 ^ 0, one can find a primal feasible dictionary by applying any ’Phase 
1 ’ algorithm that is available for the usual simplex method, see [3^ . 

The first initialization method finds a weight vector G int (7^ such that (Pi(re^)) has 
an optimal solution. Then the optimal dictionary for (Pi{w^)) is used to construct the initial 
dictionary D^. There are different ways to choose the weight vector w^. The second method 
of initialization can be thought of as a Phase 1 algorithm. It finds an initial dictionary as 
long as there exists one. 

4.8.1 Finding and constructing 

The first way to initialize the algorithm requires finding some G int such that (Pi (tu°)) 
has an optimal solution. It is clear that if the problem is known to be bounded, then any 
G int works. However, it is a nontrivial procedure in general. In the following we give 
two different methods to find such w^. The first method can also determine if the problem 
has no solution. 
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a. The first approach is to extend the idea presented in [15] to any solid polyhedral pointed 
ordering cone C. Accordingly, finding involves solving the following linear program: 


minimize 


subject to 


h^u 

\T 


(Po) 


u — Pw > 0, 

Y^{w — c) > 0, 

> 0 , 

where c G int(7^, and the columns of Y are the generating vectors of C. Under the 
assumption 5 > 0, it is easy to show that © has a maximizer if and only if ( |Pq[ ) has an 
optimal solution. Note that (Pq) is bounded. If (Pq) is infeasible, then we conclude that 


the lower image has no vertex and (|^ has no solution. In case it has an optimal solution 
(?i*,re*), then one can take G intU^, and solve (Pi(re^)) optimally. For the 

randomly generated examples of Section we have used this method. 


b. Using the idea provided in [27|, it might be possible to initialize the algorithm without 
even solving a linear program. By the structure of a particular problem, one may start 
with a dictionary which is trivially optimal for some weight . In this case, one can start 
with this choice of and get the initial dictionary even without solving an LP. An 
example is provided in Section see Remark |5.3[ 


In order to initialize the algorithm, can be used to construct the initial dictionary . 
Without loss of generality assume that — 1, indeed one can always normalize since 

c G intU implies > 0. Then, clearly G riVF. Let and be the set of basic and 

nonbasic variables corresponding to the optimal dictionary £)* of (Pi(ie^)). If one considers 
the dictionary for (P^) with the basic variables and nonbasic variables A/”^, the objective 
function of will be different from T)* as it depends on the parameter A. However, the matri¬ 
ces and hence the corresponding basic solution are the same in both dictionaries. 

We consider as the initial dictionary for the parametrized problem (Pa)- Note that is 
a nonsingular matrix as it corresponds to dictionary T)*. Moreover, since is an optimal 

dictionary for (Pi(re^)), is clearly primal feasible for (Pa) for any A G Furthermore, 

the optimality region of satisfies A^° H int A ^ 0 as A^ := [rej,..., G A^° H int A. 

Thus, is also dual feasible for (Pa) for A G A^°, and x^ is a maximizer to 0- 


4.8.2 Perturbation method 

The second method of initialization works similar to the idea presented for Algorithm [T] itself. 

Assuming that b > 0, problem (Pa) is perturbed by an additional parameter /i G M as 
follows: 


maximize (re(A)^P^ — /il^)x (Pa,/^) 

subject to Ax < 6, 

X > 0, 

where 1 is the vector of ones. After introducing the slack variables, consider the dictionary 
with basic variables ..., x^j^n ^md with nonbasic variables xi,..., This dictionary 
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is primal feasible as 6 > 0. Moreover, it is dual feasible if Pw{X) — /il < 0. We introduce the 
optimality region of this dictionary as 

:= {(A, m) e A X M+l Pw{X) - < 0}. 

Note that is not empty as n can take snfficiently large values. 

The aim of the perturbation method is to find an optimality region M such that 


Mnri(A X {0}) ^0, 


( 12 ) 


where A x {0} := {(A,0)| A G A}. If the current dictionary satisfies (12), then it can be 
taken as an initial dictionary for Algorithm [l] after deleting the parameter /i. Otherwise, 
the defining inequalities of the optimality region are found. Clearly, they correspond to the 
entering variables of the current dictionary. The search for an initial dictionary continues 
similar to the original algorithm. Note that if there does not exist a leaving variable for an 
entering variable, (Pa,/^) is found to be unbounded for some set of parameters. The algorithm 


continues until we obtain a dictionary for which the optimality region satisfies ( 12 ) or until 
we cover the the parameter set A x by the optimality regions and by the regions that 
are known to yield unbounded problems. At termination, if there exist no dictionary that 
satisfies ( 12 ), then we conclude that there is no solution to problem 0 . Otherwise, we 
initialize the algorithm with . See Example |5.l[ Remark 5.4 


5 Illustrative examples 

We provide some examples and numerical results in this section. The first example illustrates 
how the different methods of initialization and the algorithm work. The second example 
shows that Algorithm can find a solution even though the lower image does not have any 
vertices. 


Example 5.1. Consider the following problem 

maximize (xi, X 2 — X 3 , ^ 3 )^ with respect to 
subject io xi + X 2 <5 

xi + 2x2 — X2, 

X1,X2,X3 > 0. 


Let c = (1,1,1)^ G intM^. Clearly, we have A = {A G 
Let us illustrate the different initialization methods. 


Ai + A 2 ^ 1, A^ > 0, i — 1,2}. 


Remark 5.2. (Initializing by solving (Pq)? see Section |4.8.1 a.) A solution of (Pq) 
is found as re* = (1,1,1)^. Then, we take as the initial weight vector. 

x^ — (5,0,0)^ is an optimal solution found for Pi{w^). The indices of the basic variables of 
the corresponding optimal dictionary are = {I 5 5}. We form dictionary of problem (Pa) 
with basic variables B^: 


C = 5Ai 
xi = 5 

X 5 = 4 


-A1X4 
—X 4 
+X 4 


-(Ai — A 2 )x 2 

-X2 

-X2 


-(Ai + 2 A 2 — 1)X3 


+X 3 
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Remark 5.3. (Initializing using the structure of the problem, see Section 4.8.1| b.) 

The structure of Example |5 .1 1 allows us to initialize without solving a linear program. Consider 
= ( 1 , 0 , 0 )^. As the objective of Pi(u;^) is to maximize xi and the most restraining 
constraint is +X 2 <5 together with Xi > 0, x = (5, 0,0)^ is an optimal solution of Pi(u;^). 
The corresponding slack variables are X 4 = 0 and X 5 = 4. Note that this corresponds to 
the dictionary with basic variables {1,5} and nonbasic variables {2,3,4}, which yields the 
same initial dictionary as above. Note that one needs to be careful as ^ int but 
G bdC^. In order to ensure that the corresponding initial solution is a maximizer and not 
only a weak maximizer, one needs to check the optimality region of the initial dictionary. If 
the optimality region has a nonempty intersection with int A, which is the case for then 
the corresponding basic solution is a maximizer. In general, if one can find w G int such 
that ( |Pi(u;)[ ) has a trivial optimal solution, then the last step is clearly unnecessary. 


Remark 5.4. (Initializing using the perturbation method, see Section 4.8.2[ ) The 

starting dictionary for (Pa,/^) of tho perturbation method is given as 

C = -(/i-Ai)xi -(/i-A2)x2 -(/i + Ai + 2 A 2 - 1)^3 

X 4 =5 —Xi —X2 

X5 =9 —xi —2x2 +X3 

This dictionary is optimal for = {(A,/i) G AxM| fi — \i > 0, /i —A 2 > 0, /i + Ai+ 2 A 2 > 1}. 


Clearly, does not satisfy ( 12 ). The defining halfspaces for correspond to the nonbasic 
variables xi^X 2 and X 3 . If xi enters, then the leaving variable is X 4 and the next dictionary has 
the optimality region — {(A,/i) G AxM| — /i + Ai > 0, Ai — A 2 > 0, /i + Ai + 2A2 > 1} which 
satisfies (12). Then, by deleting /i the initial dictionary is found to be as above. Different 


choices of entering variables in the first iteration might yield different initial dictionaries. 


Consider the initial dictionary D^. Clearly, = {A G M^| Ai > 0}, I 2 = {A G 
M^l Ai — A 2 > 0}, and = {A G M^| Ai + 2 A 2 > 1}. The defining halfspaces for A^° H A 
correspond to the nonbasic variables X 2 , X 3 , thus we have = {2, 3}. 

The iteration starts with the only boundary dictionary D^. If X 2 is the entering variable, 
X 5 is picked as the leaving variable. The next dictionary, has basic variables = {I 52 }, 

the basic solution x^ = (1,4,0)^, and a parameter region A^ = {A G M^| 2Ai — A 2 > 
0, —Ai + A 2 > 0, 2Ai + A 2 > 1}. The halfspaces corresponding to the nonbasic variables Xj, 
j G — {3,4,5} are defining for the optimality region A^^ H A. Moreover, — {(5,2)} 
is an explored pivot for D^. 

From dictionary for entering variable xs there is no leaving variable according to the 
minimum ratio rule We conclude that problem (Pa) is unbounded for A G such that 
Ai + 2 A 2 < 1. Note that 

^11 0 


B-^N = 


-1 1 -1 


and the third column corresponds to the entering variable X 3 . 


Thus, x^o = (xJ',X 5 )^ = 


h 

(0,-1)^ and x^o = = 63 = (0,0,1)^. Thus, we add x^ = (x^^x^^x^) = 

(0, 0,1)^ to the set see Algorithm]^ line 12 . Also, by Proposition [4^ = (0, —1,1)^ 

is an extreme direction of the lower image V. After the first iteration, we have VS — {D^}, 
and BD = {D^- 
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For the second iteration, we consider G BD. There are three possible pivots with 
entering variables Xj, j G = {3,4,5}. For X 5 , X 2 is found as the leaving variable. As 
(5,2) G , the pivot is already explored and not necessary. For X 3 , the leaving variable 
is found as xi. The resulting dictionary has basic variables — {2,3}, basic solution 
x^ = (0,5,1)^, the optimality region H A = {A G M^| 2Ai + A 2 < 1, 2Ai + 3 A 2 < 
2, —Ai — 2 A 2 < —1}, and the indices of the entering variables — {1,4,5}. Moreover, we 
write — {(1,3)}. 

We continue the second iteration by checking the entering variable X 4 from . The 
leaving variable is found as xi. This pivot yields a new dictionary with basic variables 

— {2,4} and basic solution x^ — (0,4.5,0)^. The indices of the entering variables are 
found as — {1,3}. Also, we have E^^ = {(1,4)}. At the end of the second iteration we 
have VS = and BD = 

Consider G BD for the third iteration. For xi, the leaving variable is x^ and we 
obtain dictionary D^ which is already visited. For X 4 , the leaving variable is X 3 . We obtain 
the boundary dictionary D^ and update the explored pivots for it as E^^ — {(1,4), (3,4)}. 
Finally, for entering variable X 5 there is no leaving variable and one finds the same x^ that is 
already found at the first iteration. At the end of third iteration, VS = {D^^ D^^D^}^ BD = 

For the next iteration, D^ is considered. The pivots for both entering variables yield 
already explored ones. At the end of this iteration there are no more boundary dictionaries 
and the algorithm terminates with VS = {D^ ^D^ ^D^ ^D^}. Figure shows the optimality 
regions after the four iterations. The color blue indicates that the corresponding dictionary 
is visited, yellow stands for boundary dictionaries and the gray region corresponds to the set 
of parameters for which problem (P^) is unbounded. 





Figure 1: Optimality regions after the first four iterations of Example 


5.1 


The solution to the problem is (A, X^) where X = {(5,0, 0)^Al, 4, 0)^, (0, 5,1)^, (0,4.5, 0)^}, 
and X^ = {(0, 0,1)^} . The lower image can be seen in Figure^ 

Example 5.5. Consider the following example. 


maximize {xi — X 2 , X 3 — ^ 4 )^ with respect to 
subject to — X 2 + X 3 — X 4 < 1 

Xl,X2,X3,X4 > 0. 


Let c = (l,!)^ ^ intM^. Clearly, we have A = [0,1] C M. Using the method de- 
a. we find = {\^\) as the initial scalarization parameter. Then, 


scribed in Section 


4.8.1 
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Figure 2: Lower image V of Example 


5.1 


^0 _ 


= (1,0,0,0)^ is an optimal solution to Pi(u;^) and the index set of the basic vari¬ 
ables of is found as = {1}. Algorithm [b terminates after two iterations and yields 
X — {(1,0, 0,0)^, (0,0,1,0)^}, X^ — {(1, 0, 0,1) , (0,1,1,0)^}. The lower image can be seen 
in Figure Note that as it is possible to generate the lower image only with one point 
maximizer, the second one is redundant, see Remark |4.9|b. 


^ 1 
I 

CO 

0 



0 1 


Xi - X2 


2 


Figure 3: Lower image V of Example 


5.5 


6 Comparison of different simplex algorithms for LVOP 

As briefly mentioned in Section there are different simplex algorithms to solve LVOPs. 
Among them, the Evans-Steuer algorithm m works very similar to the algorithm provided 
here. It moves from one dictionary to another where each dictionary gives a point maximizer. 
Moreover, it finds ’unbounded efficient edges’, which correspond to the direction maximizers. 
Even though the two algorithms work in a similar way, they have some differences that 
affect the efficiency of the algorithms significantly. The main difference is that the Evans- 
Steuer algorithm finds the set of all maximizers whereas Algorithm finds only a subset of 
maximizers, which generates a solution in the sense of Lohne m and allows to generate the 
set of all maximal elements of the image of the feasible set. In general, the Evans-Steuer 
algorithm visits more dictionaries than Algorithm especially if the problem is degenerate. 

First of all, in each iteration of Algorithmj^ for each entering variable Xj, only one leaving 
variable is picked among the set of all possible leaving variables, see line 12. Differently, 
the Evans-Steuer algorithm performs pivots xj ^ xi for all possible leaving variables, i G 
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argmin^^^^ ^B-^N)ij>o ‘ problem is degenerate, this procedure leads the Evans- 

Steuer algorithm to visit many more dictionaries than Algorithm does. In general, these 
additionally visited dictionaries yield maximizers that are already found. In [H[2], it has been 
shown that using the lexicographic rule to choose the leaving variables would be sufficient to 
cover all the efficient basic solutions. For the numerical tests that we run, see Section we 
have modified the Evans-Steuer algorithm such that it uses the lexicographic rule. 

Another difference between the two simplex algorithms is at the step where the entering 
variables are selected. In Algorithm the entering variables are the ones which correspond 
to the defining inequalities of the current optimality region. Different methods to find the 
entering variables are provided in Section The method that is employed for the numerical 
tests of Section involves solving sequential LP’s with q — 1 variables and at most n -\- k 
inequality constraints, where k is the number of generating vectors of the ordering cone. Note 
that the number of constraints are decreasing in each LP as one solves them successively. For 
each dictionary, the total number of LPs to solve is at most n in each iteration. 

The Evans-Steuer algorithm finds a larger set of entering variables, namely ’efficient non- 
basic variables’ for each dictionary. In order to find this set, it solves n LPs with n + q + 1 
variables, q equality and n -\- q -\- 1 non-negativity constraints. More specifically, for each 
nonbasic variable j G A/* it solves 

maximize l^v 

subject to Zj^y — SZj^e^ — 'y = 0, 

y,6,v> 0, 

where y G G M,'u G M^. Only if this program has an optimal solution 0, then xj is 

an efficient nonbasic variable. This procedure is clearly costlier than the one employed in 
Algorithm]^ In [T7], this idea is improved so that it is possible to complete the procedure 
by solving fewer LPs of the same structure. Further improvements are done also in [9]. 
Moreover, in mm a different method is applied in order to find the efficient nonbasic variables. 
Accordingly, one needs to solve n LPs with 2 q variables, n equality and 2 q nonnegativity 
constraints. Clearly, this method is more efficient than the one used for the Evans-Steuer 
algorithm. However, the general idea of finding the efficient nonbasic variables clearly yields 
visiting more redundant dictionaries than Algorithmwould visit. Some of these additionally 
visited dictionaries yield different maximizers that map into already found maximal elements 
in the objective space, see Example |6.1t while some of them yield non-vertex maximal elements 
in the objective space, see Example |6.2[ 

Example 6.1. Consider the following simple example taken from [28], in which it has been 
used to illustrate the Evans-Steuer algorithm. 

maximize {3xi -h X 2 , Sxi — ^ 2 )^ with respect to 
subject to xi + X 2 < 4: 

xi — X2 < 4 

X 3 <4 

Xi,X2,X3 > 0 . 

If one uses Algorithm the solution is provided right after the initialization. The initial 
set of basic variables can be found as = {1,5,6}, and the basic solution corresponding 
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to the initial dictionary is = (4,0,0)^. One can easily check that is optimal for all 
A G A. Thus, Algorithm stops and returns the single maximizer. On the other hand, it is 
shown in [28] that the Evans-Steuer algorithm terminates only after performing another pivot 
to obtain a new maximizer x^ — (4,0,4)^. This is because, from the dictionary with basic 
variables — {1,5,6} it finds x^ as an efficient nonbasic variable and performs one more 
pivot with entering variable x^. Clearly the image of x^ is again the same vertex (4,4)^ in 
the image space. Thus, in order to generate a solution in the sense of Definition |3.1[ the last 
iteration is unnecessary. 

Example 6.2. Consider the following example. 

maximize {—xi — X 3 , —X 2 — 2 x 3 )^ with respect to 
subject to — xi — X 2 — 3 x 3 < — 1 

Xi,X2,X3 > 0. 

First, we solve the example by Algorithm Clearly, A = [0,1] C M. We find an initial 
dictionary with B^ = {1}, which yields the maximizer x^ = (1,0,0)^. One can easily 
see that index set of the defining inequalities of the optimality region can be chosen either 
as = { 2 } or = {3}. Note that Algorithm 0 picks one of them and continues with 
it. In this example we get = {2}, perform the pivot X 2 ^ xi to get with B^ = {2} 
and x^ — (0,1,0)^. From there are two choices of sets of entering variables and we 
set = {1}. As the pivot xi CG X 2 is already explored, the algorithm terminates with 
X = and X^ — 0. 

When one solves the same problem by the Evans-Steuer algorithm, from both X 2 and 
X 3 are found as entering variables. When X 3 enters from one finds a new maximizer 
x^ — (0,0, Note that this yields a nonvertex maximal element on the lower image, see 
Figure]^ 



-^2 


-1 0 

-Xi - Xs 


Figure 4: Lower image V of Example 6.2 


Remark 6.3. Note that if the problem is primal nondegenerate, then for a given entering 
variable of a given dictionary, both Algorithm [l] and the Evans-Steuer algorithm find the 
unique leaving variable. If in addition, every efficient nonbasic variable of a given dictionary 
corresponds to a defining inequality of its optimality region, then the entering variables from 
that dictionary would be the same for both algorithms. Indeed, the different type of redun¬ 


dancies that are explained in Remark 4.9 are mostly observed if there is a primal degeneracy 
or if there are efficient nonbasic variables which corresponds to redundant inequalities of the 
optimality region. Hence, it wouldn’t be wrong to state that for ’nondegenerate’ problems, the 
Evans-Steuer algorithm and Algorithm follow similar paths. But for degenerate problems 
their performance will be quite different. 
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Apart from the Evans-Steuer algorithm Ehrgott, Puerto and Rodriguez-Chia m devel¬ 
oped a primal-dual simplex algorithm to solve LVOPs. The algorithm finds a partition (A^) 
of A. It is similar to Algorithm in the sense that for each parameter set A^, it provides an 
optimal solution to the problems (P^) for all A E A^. The difference between the two algo¬ 
rithms is in the method of finding the partition. The algorithm in m starts with a (coarse) 
partition of the set A. In each iteration it finds a finer partition until no more improvements 
can be done. In contrast to the algorithm proposed here, the algorithm in [13] requires solving 
in each iteration an LP with n -\- m variables and I constraints where m < I < m -\- which 
clearly makes the algorithm computationally much more costly. In addition to solving one 
’large’ LP, it involves a procedure which is similar to finding the defining inequalities of a 
region given by a set of inequalities. Also, different from Algorithm it finds only a set of 
weak maximizers so that as a last step one needs to perform a vertex enumeration in order 
to obtain a solution consisting of maximizers only. Finally, the algorithm provided in m 
can deal with unbounded problems only if the set A^ is provided, which requires a Phase 1 
procedure. 

7 Numerical results 

In this section we provide numerical results to study the efficiency of Algorithm We 
generate random problems, solve them with different algorithms and compare the solutions 
and the CPU times. Algorithm is implemented in MATLAB. We also use a MATLAB 
implementation of Benson’s algorithm, namely bensolve 1.2 i2ni. The current version of 
bensolve 1.2 solves two linear programs in each iteration. However, we employ an improved 
version which solves only one linear program in each iteration, see [KlEaEa]. For the Evans- 
Steuer algorithm, instead of using ADBASE lai, we implement the algorithm in MATLAB. 
This way, we can test the algorithms with the same machinery. This gives the opportunity to 
compare the CPU times. For each algorithm the linear programs are solved using the GLPK 
solver, see [25] . 

The first set of problems are randomly generated with no special structure. That is to 
say, these problems are not designed to be degenerate. In particular, each element of the 
matrices A and P and the vector b is sampled independently, the elements of A and P from 
a normal distribution with mean 0 and variance 100, and the elements of b from a uniform 
distribution over [0,10]. As 5 > 0, we did not employ a Phase 1 algorithm to find a primal 
feasible initial dictionary. Table shows the numerical results for the randomly generated 
problems with three objectives. We fix different numbers of variables (n) and constraints 
(m) and generate 100 problems for each size. We measure the average time that Algorithm 
(avg A), bensolve 1.2. (avg B) and the Evans-Steuer algorithm (avg E) take to solve the 
problems. Moreover, we report the minimum (min A, min B, min E) and maximum (max A, 
max B, max E) running times for each algorithm among those 100 problems. The number of 
unbounded problems that are found among the 100 problems is denoted by #u. 

Next, we randomly generate problems with four objectives and with different numbers of 
variables (n) and constraints (m). For each size we generate four problems. Table shows 
the numerical results, where \A\ and |A^| are the number of elements of the set of point and 
direction maximizers, respectively. For each problem the time for Algorithm for bensolve 
1.2 and for the Evans-Steuer algorithm to terminate are shown by ’time A’, ’time B’ and 
’time E’, respectively. 


23 


Table 1: Run time statistics for randomly generated problems where q — ?>. For the first row 
n — 20, m = 40; for the second row, n = 30, m — 30; for the last row n = 40, m = 20. 


min A 

min B 

min E 

max A 

max B 

max E 

avg A 

avg B 

avg E 

#u 

0.34 

0.20 

0.27 

6.02 

162.53 

5.77 

1.72 

15.63 

1.63 

0 

0.08 

0.08 

0.09 

9.36 

257.98 

8.61 

3.15 

32.41 

2.90 

8 

0.05 

0.03 

0.08 

13.52 

418.44 

11.81 

3.33 

23.92 

2.96 

38 


For these particular examples all algorithms find the same solution (T, X^). As no struc¬ 
ture is imposed on these problems, the probability that these problems are nondegenerate is 
very high. This explains finding the same solution by all of the algorithms. As seen from 
the Tables and the CPU times of the Evans-Steuer algorithm are very close to the CPU 
times of Algorithm which is expected as explained in Remark |6.3[ 

Table 2: Computational results for randomly generated problems 


q 

n 

m 

|;e| 

\X^\ 

time A 

time B 

time E 

4 

30 

50 

267 

0 

3.91 

64.31 

3.61 

4 

30 

50 

437 

0 

6.95 

263.39 

7.06 

4 

30 

50 

877 

0 

15.73 

1866.1 

17.01 

4 

30 

50 

2450 

0 

74.98 

33507 

73.69 

4 

40 

40 

814 

0 

20.41 

1978.3 

18.56 

4 

40 

40 

1468 

81 

42.39 

11785 

38.20 

4 

40 

40 

2740 

0 

105.45 

64302 

97.69 

4 

40 

40 

2871 

324 

121.16 

82142 

112.11 

4 

50 

30 

399 

21 

10.53 

233.11 

9.23 

4 

50 

30 

424 

0 

11.22 

294.17 

9.92 

4 

50 

30 

920 

224 

28.08 

3434.1 

24.05 

4 

50 

30 

1603 

176 

55.97 

14550 

49.86 


As seen from Tables and the parametric simplex algorithm works more efficiently 
than bensolve 1.2 for the randomly generated problems. The main reason for the difference 
in the performances is that in each iteration, Benson’s algorithm solves an LP that is in the 
same size of the original problem and also a vertex enumeration problem. Note that solving a 
vertex enumeration problem from scratch in each iteration is a costly procedure. In um , an 
online vertex enumeration method has been proposed and this would increase the efficiency 
of Benson’s algorithm. 

Note that these randomly generated problems have no special structure and thus there 
is a high probability that these problems are nondegenerate. However, in general, Benson- 
type objective space algorithms are expected to be more efficient whenever the problem is 
degenerate. The main reason is that these algorithms do not need to deal with the different 
efficient solutions which map into the same point in the objective space. This, indeed, is one 
of the main motivation of Benson’s algorithm for linear multiobjective optimization problems, 
see [5]. 
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In order to see the efficiency of our algorithm for degenerate problems, we generate random 
problems which are designed to be degenerate. In the following examples this is done by gen¬ 
erating a nonnegative b vector with many zero components and choosing objective functions 
with the potential to create optimality regions with empty interior within A. In particular, 
for the three-objective examples, we generate the first objective function randomly, take the 
second one to be the negative of the first objective function, and let the third objective con¬ 
sist of only one nonzero entry. For the four-objective examples, the first three objectives are 
created as described above and the fourth one is generated randomly in a way that at least 
half of its components are zero. The number of nonzero elements in b and in the last objective 
function of the four-objective problems are sampled independently from uniform distributions 
over the integers in the intervals [0, [^J] and [0, [|J], respectively. Each element of the first 
column and each possibly nonzero element of the third and the fourth column of P as well 
as each element of A and each possibly nonzero element of b is sampled independently, in the 
same way as for the nondegenerate problems. 

First, we consider three objective functions where we fix different numbers of variables (n) 
and constraints (m). We generate 20 problems for each size. We measure the average time 
that Algorithmic (avg A), bensolve 1.2. (avg B) and the Evans-Steuer algorithm (avg E) take 
to solve the problems. We also report the minimum (min A, min B, min E) and maximum 
(max A, max B, max E) running times for each algorithm among those 20 problems. The 
times are measured in seconds. The results are given in Table 

Table 3: Run time statistics for randomly generated degenerate problems where q = 3. For the 
first row n = 5, m = 15; for the second row n = m = 10; and for the last row n = 15, m = 5. 


min A 

min B 

min E 

max A 

max B 

max E 

avg A 

avg B 

avg E 

0.02 

0.01 

0.03 

0.14 

0.09 

140.69 

0.07 

0.04 

8.85 

0.03 

0.02 

0.08 

1.33 

0.14 

4194.1 

0.25 

0.04 

227.02 

0.05 

0.01 

0.09 

1.47 

0.20 

1893.0 

0.05 

0.25 

190.25 


In order to give an idea how the solutions provided by the three algorithms differ for these 
degenerate problems, in Table we provide detailed results for single problems. Among the 
20 problems that are generated to obtain each row of Table we select the two problems 
with the CPU times ’max A’ and ’max E’ and provide the following for them. |T(.)| and 
denote the number of elements of the set of point and direction maximizers that are 
found by each algorithm, respectively. IUS^aI and lUS^I are the number of dictionaries that 
Algorithm [T] and the Evans-Steuer algorithm visit until termination. For each problem the 
time for Algorithm [C bensolve 1.2, and the Evans Steuer algorithm to terminate are shown 
by ’time A’, ’time B’ and ’time E’, respectively. 

Finally, we compare Algorithmic and bensolve 1.2 to get statistical results regarding their 
efficiencies for degenerate problems. Note that this test was done on a different computer 
than the previous tests. Table |C shows the numerical results for the randomly generated 
degenerate problems with g = 4 objectives, m constraints and n variables. We generate 100 
problems for each size. We measure the average time that Algorithm |C (avg A) and bensolve 
1.2. (avg B) take to solve the problems. The minimum (min A, min B) and maximum (max 
A, max B) running times for each algorithm among those 100 problems are also provided. 

Clearly, for degenerate problems bensolve 1.2 is more efficient than the simplex-type al- 
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Table 4: Computational results for single problems that require CPU times max A and max 
E among the ones that are generated for TableFor the first set of problems n = 5, m = 15; 
for the second set of problems n = m = 10 (max A and max E yielded the same problem 
here); and for the last set of problems n = 15, m = 5. 



\VSe\ 


\^b\ 

|Ae| 

1 'yh\ 


1 'yh 1 

time A 

time B 

time E 

20 

5617 

1 

1 

1 

0 

0 

0 

0.13 

0.03 

140.69 

30 

361 

3 

3 

3 

0 

0 

0 

0.14 

0.06 

6.61 

324 

22871 

1 

1 

4 

0 

0 

1 

1.33 

0.03 

4194.1 

14 

11625 

1 

1 

1 

0 

0 

0 

0.05 

0.03 

1893.0 

452 

11550 

1 

1 

1 

39 

2 

4707 

1.47 

0.05 

1598.1 


Table 5: Run time statistics for randomly generated degenerate problems. 


q 

n 

m 

min A 

min B 

max A 

max B 

avg A 

avg B 

4 

10 

30 

0.05 

0.03 

177.47 

4.07 

8.49 

0.21 

4 

20 

20 

0.21 

0.02 

973.53 

199.77 

19.89 

12.05 

4 

30 

10 

0.11 

0.03 

2710.20 

13.70 

37.68 

0.73 


gorithms considered here, namely Algorithm and the Evans-Steuer algorithm. However, 
the design of Algorithm results in a significant decrease in CPU time compared to the 
Evans-Steuer algorithm in its improved form of mm- 
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